
cd "C:\Users\Administrator\Desktop\zhouxb"

*------------------
*- 方法1：连玉君
*------------------


clear 

set obs 1000
set seed 123
gen x = rnormal()
gen v = rnormal(0.1, 1.2)
save a.dta, replace


/*
clear all

input x     v 
      0.11  0.1
      0.22  0.2
      0.33  0.3
end
save a.dta, replace 
*/
timer on 1 

preserve 
  use a.dta, clear 
  gen id = _n
  keep id x 
  save ax.dta, replace
restore

preserve 
  use a.dta, clear 
  keep v
  save av.dta, replace
restore

use ax.dta, clear 
cross using av.dta

sort id

gen xvi = x+v
gen double Norm_xvi = normal(xvi)
bysort id: egen mean_xv = mean(Norm_xvi)      // New
bysort id: egen mean_xv_w = mean(xvi*Norm_xvi)
duplicates drop id, force
list x v mean*

timer off 1

timer list 1
 
 
 
*------------------
*- 方法2：江艇 
*------------------


clear all

input x     v 
      0.11  0.1
      0.22  0.2
      0.33  0.3
end
gen ID = _n
save a.dta, replace 


/*
*-需要预先安装 moremata 命令
  *-方法1：官方安装，很慢
    ssc install moremata, replace   
  *-方法2：使用 tssc 命令安装 (推荐)
    net install tssc.pkg, from("https://tidyfriday.gitee.io/tssc/")
    tssc install moremata, replace
*/

clear
set obs 80000
set seed 123
gen x = rnormal()
gen v = rnormal(0.1, 1)
gen ID = _n // 观察值序号
save a.dta, replace


use "a.dta", clear
global N = _N   // 样本数

timer on 2       // 计时器

mata: mata clear // 清空内存中的 mata 对象
putmata x v //ID // 把变量转换为 mata 矩阵


*-修正了江艇代码的小错误，他把 x 和 v 搞反了
mata
    v  = v'
    vv = mm_expand(v, $N) 
    xv = x:+vv
//xv	
    xv2= normal(xv)
//xv2	
    mean_xv2 = mean(xv2',1)  // column means
//mean_xv2
    mean_xv2 = mean_xv2'
    xv3 = xv:*xv2
    mean_xv3 = mean(xv3',1)
    mean_xv3 = mean_xv3'
end
getmata mean_xv2 mean_xv3 // , id(ID)  // force

list in 1/3

timer off  2
timer list 2  

*-耗时测试 
* N			耗时
* 10000		 7.02 s
* 50000		43.26 s
* 100000    出错了

* 连玉君：王俊，看起来，江艇的代码不存在你说的 N --> N^2 问题
*         所以，我目前还是没想明白哪里出错了？
/*
    v  = v'               // 1xN
    vv = mm_expand(v, $N) // NxN 
*/

*- Mata Limits:  -help m1 limits-   这些内容放到第一小节的末尾
/*
								min		maximum
_________________________________________________________				
Scalars, vectors, matrices
	rows 				  		0 		2,147,483,647
	columns 			  		0 		2,147,483,647
String elements, length			0 		2,147,483,647
_________________________________________________________
Stata’s matsize plays no role in these limits.
*/



*------------------
*- 江艇改进
*------------------


clear all

input x     v 
      0.11  0.1
      0.22  0.2
      0.33  0.3
end
gen ID = _n
save a.dta, replace 

use "a.dta", clear
global N = _N   // 样本数

timer on 3       // 计时器

mata: mata clear // 清空内存中的 mata 对象

putmata x v // 把变量转换为 mata 矩阵

mata
    x  = x'
    xx = mm_expand(x, $N) 
  xx
    xv = xx:+v
  xv
    //xv2= normal(xv)
	xv2 = xv
    mean_xv2 = mean(xv2,1)
  mean_xv2
    mean_xv2 = mean_xv2'
    xv3 = xv:*xv2
    mean_xv3 = mean(xv3,1)
    mean_xv3 = mean_xv3'
end

getmata mean_xv2 mean_xv3 // force


timer list 3
